---
output:
  html_document: default
  pdf_document: default
  word_document: default
---

```{r}

##Creating Status Loss Replication Code##
##Author: Alex Yu-Ting Lin, https://sites.google.com/view/alexytlin/
##Version date: May 3, 2025
##Estimate run time: 3-4 hours


##loading libraries
library(quanteda)
library(tidytext)
library(topicmodels)
library(tidyverse)
library(ggplot2)
library(scales)
library(ldatuning)
library(stopwords)
library(writexl)
#library(easyGgplot2)
library(stm)
library(dplyr)
#library(padr)
#library(reshape)
library(gridExtra)
library(jtools)
library(stargazer)
library(ggstance)
library(sqldf)
library(broom.mixed)
library(stm)
library(survey)
library(MASS)
library(statmod)
library(mediation)
library(vtable)
library(ggpubr)
library(rstatix)
library(pwr)
library(syuzhet)
library(ggfortify)


##loading data

#main survey data
data <- read.csv("/creatingstatuslossdata.csv", stringsAsFactors= F)

#text data for Appendix 2
mfa <- read.csv("/mfa.csv", stringsAsFactors= F)


##cleaning survey data

data$treatment <- ifelse(
  (data$group == "Treatment"), 1, 0)

##creating demographic variables

data$male <- ifelse((data$gender == "2"), 1, 0)
data$yob <- as.numeric(data$yob)
data$age <- 2021 - data$yob


##creating DV variables

data$rankingscore <- (data$now + data$tenyears + data$shouldlead + data$deservelead)/4
data$respectscore <- (data$ownperception + data$otherperception + data$respect + data$morality + data$principles + data$emulation) /6 
data$competencescore <- (data$rightthing + data$better + data$domore) /3
data$statusscore <- (data$rankingscore + data$respectscore + data$competencescore) /3

```

```{r}

#########################
###MODELS FOR TABLE 2####
#########################

rankingmodel <- lm(formula = rankingscore ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)

respectmodel <- lm(formula = respectscore ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)

competencemodel <- lm(formula = competencescore ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)

#######################
###CREATING TABLE 2####
#######################

stargazer(rankingmodel, respectmodel, competencemodel, covariate.labels=c("Delegitimation (treatment)","Age",
"Male","Education","Political Leaning","Income", "Knowledge of International Affairs"), type = "latex")

```

```{r}

#############################################
#####CREATING GRAPH 2 BY BREAKING DVS########
#############################################

##RANKING QUESTIONS##
nowmodel <- lm(formula = now ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)
tenyearsmodel <- lm(formula = tenyears ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)
shouldleadmodel <- lm(formula = shouldlead ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)
deserveleadmodel <- lm(formula = deservelead ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)

##RANKING ROBUSTNESS VISUALIZATION##
rankingrobustness <- plot_summs(nowmodel, tenyearsmodel, shouldleadmodel, deserveleadmodel, scale = TRUE, ci_level = 0.95, inner_ci_level = .834, coefs = c("Delegitimation (treatment)" = "treatment"), model.names = c("US ranking now", "US ranking in ten years", "Should US lead?", "Does the US deserve to lead?"))

##RESPECT QUESTIONS##
ownperceptionmodel <- lm(formula = ownperception ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)
otherperceptionmodel <- lm(formula = otherperception ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)
respectrmodel <- lm(formula = respect ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)
moralitymodel <- lm(formula = morality ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)
principlesmodel <- lm(formula = principles ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)
emulationmodel <- lm(formula = emulation ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)

##RESPECT ROBUSTNESS VISUALIZATION##
respectrobustness <- plot_summs(ownperceptionmodel, otherperceptionmodel, respectrmodel, moralitymodel, principlesmodel, emulationmodel, scale = TRUE, ci_level = 0.95, inner_ci_level = .834, coefs = c("Delegitimation (treatment)" = "treatment"), model.names = c("Own perception of US", "Others' perception of US", "Do people respect the US?", "Does the US act on morality??", "Does the US act on principles? ", "Is the US a model for emulation?"))

##COMPETENCE QUESTIONS

rightthingmodel <- lm(formula = rightthing ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)
bettermodel <- lm(formula = better ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)
domoremodel <- lm(formula = domore ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)

##COMPETENCE ROBUSTNESS VISUALIZATION##

competencerobustness <- plot_summs(rightthingmodel, bettermodel, domoremodel, deserveleadmodel, scale = TRUE, inner_ci_level = .834, coefs = c("Delegitimation (treatment)" = "treatment"), model.names = c("Do the right thing?", "Does the US make the world better?", "Should the US do more?"))


##CALLING GRAPH 2 IN MAIN MANUSCRIPT

graph2 <- ggarrange(rankingrobustness, respectrobustness, competencerobustness,
ncol = 1, nrow = 3)

graph2

##SAVING GRAPH 2

```

```{r}

#############
##MEDIATION##
#############

#1. Military exercise mediation 

#delegitimation -> exercises
mildirectmodel <- lm(formula = exercise ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)

#statusmodel <- lm(formula = statusscore ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)

#delegitimation -> status 
med.status <- lm(statusscore ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)

#status -> exercises 
mil.status <- lm(exercise ~ statusscore + treatment + age + male + education + leaning + income + internationalaffairs, data = data)

#mediation effects 
milmediation <- mediate(med.status, mil.status, treat = "treatment", mediator = "statusscore", robustSE = F , boot = TRUE, sims = 2000, conf.level = .95)

#call Figure 3 results
summary(milmediation)

```

```{r}

#2. TPP mediation

#delegitimation -> tpp
tppdirectmodel <- lm(tpp ~ treatment + age + male + education + leaning + income + internationalaffairs, data = data)

#status -> tpp 
tpp.status <- lm(tpp ~ statusscore + treatment + age + male + education + leaning + income + internationalaffairs, data = data)

#note: med.status is the same model for delegitimation -> status; med.status is provided above#

#call mediation on tpp 
tppmediation <- mediate(med.status, tpp.status, treat = "treatment", mediator = "statusscore", robustSE = F , boot = TRUE, sims = 2000, conf.level = .95)

##Call Figure 4
summary(tppmediation)

```

```{r}

#################################
###END OF MAIN MANUSCRIPT CODE###
#################################

```

```{r}

#######################
#######APPENDIX########
#######################


########################
#APPENDIX #1: DESCRIPTIVE STATISTICS#
########################

descriptivestatistics <- stargazer(data)

```
```{r}

######################################################
####APPENDIX 2: TEXT ANALYSIS FOR GENERALIZABILITY####
######################################################

##cleaning data: create variable, delete duplicates, and define corpus

#duplicated(mfa)
mfa <- mfa[!duplicated(mfa$answer), ]
##percentile
#quantile(mfa$a_sentiment, c(0, .1, .2, .3, .4, .5), na.rm = FALSE)
mfacriticism <- mfa[mfa$a_sentiment <= 2.0079,]
text <- (mfacriticism$answer)
corpus <- corpus(mfacriticism$answer)


##cleaning text

corpus <- tokens(corpus, remove_punct = TRUE, remove_numbers = TRUE)
stopwords <- c("china", stopwords("english"))
corpus <- tokens_remove(corpus, pattern = stopwords)

mfadfm <- dfm(corpus)
#mfadfm <- dfm_trim(mfadfm, min_docfreq=100) %>%
   #dfm_subset(ntoken(mfadfm) > 0)
              
topfeatures(mfadfm, 2000)


##LDA topic modeling
lda <- LDA(mfadfm, k = 30, method = "Gibbs", control = list(verbose=25L, seed = 123, burnin = 100, iter = 3000))

terms <- get_terms(lda, 40)
terms

```

```{r}
topics <- get_topics(lda, 1) ##what's highest topic of a given document? 


mfacriticism$prob_topic <- lda@gamma[,1:30]
mfacriticism$pred_topic <- topics
topicnumbers <- table(topics)

#calculating # of documents by topic, the #s used in Appendix 2 tables
topicnumbers 

```

```{r}

##Selecting topics

##Determine optimal number of topics
library(tm)
corpus <- Corpus(VectorSource(mfacriticism$answer))
testcorpus <- TermDocumentMatrix(corpus)
dtm <- testcorpus

##Method #1 from LDA Tuning

result <- FindTopicsNumber(
  dtm,
  topics = seq(from = 2, to = 60, by = 5),
  metrics = c("Griffiths2004", "CaoJuan2009", "Arun2010", "Deveaud2014"),
  method = "Gibbs",
  control = list(seed = 77),
  verbose = TRUE
)

topicnumbergraph1 <- FindTopicsNumber_plot(result)

```

```{r}
##Method #2: perplexity test

n_topics <- c(2, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50)
compare <- n_topics %>%
  map(LDA, x = dtm, control = list(seed = 123))
data_frame(k = n_topics,
           perplex = map_dbl(compare, perplexity)) %>%
  ggplot(aes(k, perplex)) +
  geom_point() +
  geom_line() +
  labs(x = "Number of topics",
       y = "Perplexity")

```

```{r}

##Text analysis validation by year

#covid validation, which happened in 2020
covid <- mfacriticism[which(mfacriticism$pred_topic==20),]
covidbyyear <- count(covid, year)

#insert 0 then plot 

covidbyyear <- rbind(covidbyyear, data.frame("year"="2002", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2003", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2004", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2005", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2006", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2007", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2008", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2009", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2010", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2011", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2012", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2013", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2014", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2015", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2016", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2017", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2018", "n"=0))
covidbyyear <- rbind(covidbyyear, data.frame("year"="2019", "n"=0))

##plot COVID-19 graph

covidgraph <- ggplot(data=covidbyyear, aes(x= year, y=n)) + geom_bar(position="dodge", stat="identity", width = 0.3, color = "black") + ylab ("Frequency") + xlab ("") + ggtitle("Coverage of COVID-19") + theme_classic()

covidgraph <- covidgraph + ylab ("Frequency") + xlab ("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1), axis.title.x = element_text(hjust= 0.5), axis.title.y = element_text(vjust= 0.5))  

## call covid validation graph
covidgraph

#hk protest validation, Umbrella 2014, then 2019 over national security bill

hkprotest <- mfacriticism[which(mfacriticism$pred_topic==22),]
hkprotestbyyear <- count(hkprotest, year)

hkprotestbyyear <- rbind(hkprotestbyyear, data.frame("year"="2004", "n"=0))
hkprotestbyyear <- rbind(hkprotestbyyear, data.frame("year"="2005", "n"=0))
hkprotestbyyear <- rbind(hkprotestbyyear, data.frame("year"="2006", "n"=0))
hkprotestbyyear <- rbind(hkprotestbyyear, data.frame("year"="2008", "n"=0))
hkprotestbyyear <- rbind(hkprotestbyyear, data.frame("year"="2012", "n"=0))
hkprotestbyyear <- rbind(hkprotestbyyear, data.frame("year"="2013", "n"=0))


hkgraph <- ggplot(data=hkprotestbyyear, aes(x= year, y=n)) + geom_bar(position="dodge", stat="identity", width = 0.3, color = "black") + ylab ("Frequency") + xlab ("") + ggtitle("Coverage of Hong Kong Protest") + theme_classic()

hkgraph <- hkgraph + ylab ("Frequency") + xlab ("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1), axis.title.x = element_text(hjust= 0.5), axis.title.y = element_text(vjust= 0.5))  

#call hk validation graph

hkgraph

#sensitivity analysis 

mfacriticism$sensitivity <- ifelse(
  (
    mfacriticism$pred_topic == "11" |
      mfacriticism$pred_topic == "13" |
      mfacriticism$pred_topic == "18" |
      mfacriticism$pred_topic == "20" |
      mfacriticism$pred_topic == "21" |
      mfacriticism$pred_topic == "23" |
      mfacriticism$pred_topic == "9" |
      mfacriticism$pred_topic == "25" |
      mfacriticism$pred_topic == "28"
  ),
  1,
  ifelse((
    mfacriticism$pred_topic == "6" |
      mfacriticism$pred_topic == "15" |
      mfacriticism$pred_topic == "19" |
      mfacriticism$pred_topic == "22" |
      mfacriticism$pred_topic == "26"
  ),
  2,
  3
  )
)

quantile(mfa$a_sentiment, c(0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3), na.rm = FALSE)

mfacriticism$sensitivity <- as.numeric(mfacriticism$sensitivity)
count(mfacriticism, sensitivity)

criticism5 <- mfacriticism[which(mfacriticism$a_sentiment <= 0.57165),]
criticism10 <- mfacriticism[which(mfacriticism$a_sentiment <= 0.94930),]
criticism15 <- mfacriticism[which(mfacriticism$a_sentiment <= 1.26100),]
criticism20 <- mfacriticism[which(mfacriticism$a_sentiment <= 1.52160),]
criticism25 <- mfacriticism[which(mfacriticism$a_sentiment <= 1.76500),]

count(criticism5, sensitivity)
count(criticism10, sensitivity)
count(criticism15, sensitivity)
count(criticism20, sensitivity)
count(criticism25, sensitivity)

```



```{r}


#################
####APPENDIX #3##
#################

#################
##BALANCE TESTS##
#################

sumtable(data, group = 'treatment', group.test = T, out = 'latex')

##########
##T-TESTS#
##########

##T-test
tranking <- t.test(rankingscore~treatment, data=data)
trespect <- t.test(respectscore~treatment, data=data)
tcompetence <- t.test(competencescore~treatment, data=data)

##ranking t-test plot##

rankingplot <- ggboxplot(data, x = "treatment", y = "rankingscore",
          ylab = "Ranking", xlab = "") +
          ylim(0,7.5)

rankingplot <- rankingplot + scale_x_discrete(labels=c("0" = "Control", "1" = "Treatment"))
rankingplot <- rankingplot + stat_compare_means(method = "t.test")


##respect t-test plot##

respectplot <- ggboxplot(data, x = "treatment", y = "respectscore",
          ylab = "Respect", xlab = "") +
          ylim(0,7.5)

respectplot <- respectplot + scale_x_discrete(labels=c("0" = "Control", "1" = "Treatment"))
respectplot <- respectplot + stat_compare_means(method = "t.test")


##competence t.test plot## 

competenceplot <- ggboxplot(data, x = "treatment", y = "competencescore",
          ylab = "Competence", xlab = "") +
          ylim(0,7.5)

competenceplot <- competenceplot + scale_x_discrete(labels=c("0" = "Control", "1" = "Treatment"))
competenceplot <- competenceplot + stat_compare_means(method = "t.test") + geom_boxplot(outlier.shape = NA)

##combining t-test plots

ttestplots <- ggarrange(rankingplot, respectplot, competenceplot,
ncol = 3, nrow = 1)

##call t-test plots
ttestplots

```


```{r}

############################
##APPENDIX 4: PCA ANALYSIS##
############################

status.pca <- prcomp(data[,c(33:45)], center = TRUE,scale. = TRUE)

pcaplot <- autoplot(status.pca, data = data, size = 0.75,
         loadings = TRUE, loadings.label = TRUE,
         loadings.label.size = 4, loadings.label.colour = 'blue',
         loadings.label.vjust = 0.75,
         loadings.label.hjust = 0.05,
         loadings.colour = 'blue')

pcaplot <- pcaplot + theme_classic()

#Call PCA plot
pcaplot

```

```{r}

#############################################
##APPENDIX 5: REGRESSIONS WITHIOUT CONTROLS## 
#############################################

#Replicating Table 2

rankingpuremodel <- lm(formula = rankingscore ~ treatment, data = data)
respectpuremodel <- lm(formula = respectscore ~ treatment, data = data)
competencepuremodel <- lm(formula = competencescore ~ treatment, data = data)
stargazer(rankingpuremodel, respectpuremodel, competencepuremodel, type = "latex")

```

```{r}
#Replicating Figure 2

#RANKING
nowpuremodel <- lm(formula = now ~ treatment, data = data)
tenyearspuremodel <- lm(formula = tenyears ~ treatment, data = data)
shouldleadpuremodel <- lm(formula = shouldlead ~ treatment, data = data)
deserveleadpuremodel <- lm(formula = deservelead ~ treatment, data = data)


##full
rankingpuremodelfull <- plot_summs(nowpuremodel, tenyearspuremodel, shouldleadpuremodel, deserveleadmodel, scale = TRUE, ci_level = 0.95, inner_ci_level = .834, coefs = c("Delegitimation (treatment)" = "treatment"), model.names = c("US ranking now", "US ranking in ten years", "Should US lead?", "Does the US deserve to lead?"))


#RESPECT
ownperceptionpuremodel <- lm(formula = ownperception ~ treatment, data = data)
otherperceptionpuremodel <- lm(formula = otherperception ~ treatment, data = data)
respectrpuremodel <- lm(formula = respect ~ treatment, data = data)
moralitypuremodel <- lm(formula = morality ~ treatment, data = data)
principlespuremodel <- lm(formula = principles ~ treatment, data = data)
emulationpuremodel <- lm(formula = emulation ~ treatment, data = data)

##full
respectrpuremodelfull <- plot_summs(ownperceptionpuremodel, otherperceptionpuremodel, respectrpuremodel, moralitypuremodel, principlespuremodel, emulationpuremodel, scale = TRUE, ci_level = 0.95, inner_ci_level = .834, coefs = c("Delegitimation (treatment)" = "treatment"), model.names = c("Own perception of US", "Others' perception of US", "Do people respect the US?", "Does the US act on morality??", "Does the US act on principles? ", "Is the US a model for emulation?"))


#COMPETENCE

rightthingpuremodel <- lm(formula = rightthing ~ treatment, data = data)
betterpuremodel <- lm(formula = better ~ treatment, data = data)
domorepuremodel <- lm(formula = domore ~ treatment, data = data)

##full

competencepuremodelfull <- plot_summs(rightthingpuremodel, betterpuremodel, domorepuremodel, deserveleadmodel, scale = TRUE, inner_ci_level = .834, coefs = c("Delegitimation (treatment)" = "treatment"), model.names = c("Do the right thing?", "Does the US make the world better?", "Should the US do more?"))


appendixpuregraph <- ggarrange(rankingpuremodelfull, respectrpuremodelfull, competencepuremodelfull,
ncol = 1, nrow = 3)

#Replicating Figure 

appendixpuregraph

```

```{r}

#Replicating Figure 3

#delegitimation -> exercises
puremildirectmodel <- lm(formula = exercise ~ treatment, data = data)


#delegitimation -> status 
puremed.status <- lm(statusscore ~ treatment, data = data)

#status -> exercises 
puremil.status <- lm(exercise ~ statusscore + treatment, data = data)

#call mediation effects 
puremilmediation <- mediate(puremed.status, puremil.status, treat = "treatment", mediator = "statusscore", robustSE = F , boot = TRUE, sims = 2000, conf.level = .95)

summary(puremilmediation, type = "latex")

#stargazer(puremed.status, puremil.status, puremildirectmodel, type = "latex")


#Replicating Figure 4 

#delegitimation -> tpp
puretppdirectmodel <- lm(formula = tpp ~ treatment, data = data)

#status -> tpp 
puretpp.status <- lm(tpp ~ statusscore + treatment, data = data)


##call mediation on tpp 
puretppmediation <- mediate(puremed.status, puretpp.status, treat = "treatment", mediator = "statusscore", robustSE = F , boot = TRUE, sims = 2000, conf.level = .95)

summary(puretppmediation, type = "latex")
#stargazer(puremed.status, puretpp.status, puretppdirectmodel, type = "latex")

```

```{r}

#############################################################
####APPPENDIX 6: FULL REGRESSION TABLE FOR MIL MEDIATION#####
#############################################################

stargazer(med.status, mil.status, mildirectmodel, type = "latex")
summary(milmediation, type = "latex")


```


```{r}

#############################################################
####APPPENDIX 7: FULL REGRESSION TABLE FOR TPP MEDIATION#####
#############################################################

stargazer(med.status, tpp.status, tppdirectmodel, type = "latex")
summary(tppmediation, type = "latex")

```


```{r}

########################################
##APPENDIX #8: MIL MEDATION SENSITIVITY#
########################################

milsensitivity <- medsens(milmediation, rho.by = 0.1, effect.type = "indirect", sims = 2000)

summary(milsensitivity)
plot(milsensitivity)

```

```{r}

#########################################
##APPENDIX #8: TPP MEDIATION SENSITIVITY#
#########################################

#tpp sensitivity 

tppsensitivity <- medsens(tppmediation, rho.by = 0.1, effect.type = "indirect", sims = 2000)
summary(tppsensitivity)
plot(tppsensitivity)

```

```{r}


###################################
###APPENDIX 9: SUBGROUP ANALYSIS###
###################################

#age, male, education, income, international affairs

subrankingmodel <- lm(formula = rankingscore ~ treatment + age + male + education + income + internationalaffairs + treatment*age + treatment*male + treatment*education + treatment*income + treatment*internationalaffairs, data = data)

subrespectmodel <- lm(formula = respectscore ~ treatment + age + male + education + income + internationalaffairs + treatment*age + treatment*male + treatment*education + treatment*income + treatment*internationalaffairs, data = data)

subcompetencemodel <- lm(formula = competencescore ~ treatment + age + male + education + income + internationalaffairs + treatment*age + treatment*male + treatment*education + treatment*income + treatment*internationalaffairs, data = data)

#Call Appendix 9 regression table
stargazer(subrankingmodel, subrespectmodel, subcompetencemodel, type = "latex")

```

```{r}

####################################
#APPENDIX 9 MIL SUBGROUP MEDIATION#
####################################

submed.status <- lm(statusscore ~ treatment*internationalaffairs + age + male + education + leaning + income, data = data)

#status -> exercises 
submil.status <- lm(exercise ~ treatment*internationalaffairs + statusscore*internationalaffairs + age + male + education + leaning + income + internationalaffairs, data = data)

submilmediation <- mediate(submed.status, submil.status, treat = "treatment", mediator = "statusscore", robustSE = F , boot = TRUE, sims = 2000, conf.level = .95)

test.modmed(submilmediation, list(internationalaffairs = 1), list(internationalaffairs = 2), sims = 2000)



```


```{r}

#####################################
#APPENDIX 9 TPP SUBGROUP MEDIATION##
#####################################

submed.status <- lm(statusscore ~ treatment*internationalaffairs + age + male + education + leaning + income, data = data)

#status -> exercises 
subtpp.status <- lm(tpp ~ treatment*internationalaffairs + statusscore*internationalaffairs + age + male + education + leaning + income + internationalaffairs, data = data)

subtppmediation <- mediate(submed.status, subtpp.status, treat = "treatment", mediator = "statusscore", robustSE = F , boot = TRUE, sims = 2000, conf.level = .95)

test.modmed(subtppmediation, list(internationalaffairs = 1), list(internationalaffairs = 2), sims = 2000)

```

